7-3 all cause mortality

Analyses of all cause mortality.

Authors

Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Load packages

library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)

Prepare datasets

all_dat <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_dat)
fas_itt_nona_dat <- fas_itt_dat |>
  filter(!is.na(D28_death))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_dat)
acs_itt_nona_dat <- acs_itt_dat |>
  filter(!is.na(D28_death))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_dat)
avs_itt_nona_dat <- avs_itt_dat |>
  filter(!is.na(D28_death))

Load models

logistic_mod <- compile_cmdstanr_mod(
  file.path("binary", "logistic"), dir = "stan")
logisitc_site <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site"), dir = "stan")
logistic_site_epoch <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")

Code

make_D28_death_table_A <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(AAssignment = factor(
      AAssignment, 
      levels = paste0("A", 0:2),
      labels = intervention_labels()$AAssignment)
    ) %>%
    group_by(AAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(AAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(AAssignment = fct_inorder(AAssignment)) %>%
    gather(key, value, -AAssignment, factor_key = T) %>%
    spread(AAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}


make_D28_death_table_C <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = paste0("C", 0:4),
      labels = intervention_labels()$CAssignment)
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}

Descriptive

Anticoagulation

Relationship anti-coagulation to outcome

make_D28_death_table_C(
    acs_itt_dat
)

n (\%)	Low dose	Intermediate dose	Low dose with aspirin	Therapeutic dose	Overall
Randomised	610	613	283	50	1556
Outcome missing	12 (2.0)	5 (0.8)	2 (0.7)	0 (0.0)	19 (1.2)
Outcome observed	598 (98.0)	608 (99.2)	281 (99.3)	50 (100.0)	1537 (98.8)
Died within 28 days	19 (3.2)	15 (2.5)	10 (3.6)	6 (12.0)	50 (3.3)

Table 1: Primary outcome by anti-coagulation intervention.

Relationship anti-coagulation to outcome

save_tex_table(make_D28_death_table_C(
      acs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-3-anticoagulation-summary")

Antiviral

Relationship anti-viral to outcome

make_D28_death_table_A(
     avs_itt_dat
)

n (\%)	Standard of care	Nafamostat	Overall
Randomised	73	82	155
Outcome missing	0 (0.0)	0 (0.0)	0 (0.0)
Outcome observed	73 (100.0)	82 (100.0)	155 (100.0)
Died within 28 days	0 (0.0)	0 (0.0)	0 (0.0)

Table 2: Primary outcome by anti-viral intervention.

Relationship anti-viral to outcome

save_tex_table(make_D28_death_table_A(
     avs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-3-antiviral-summary")

Age

Relationship age to outcome

agedat <- fas_itt_dat %>%
  dplyr::count(D28_death, AgeAtEntry) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\ndied by day 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and the primary outcome.

Sex

Relationship sex to outcome

tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(Sex, D28_death) %>%
  group_by(Sex) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Sex, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants")
p2 <- ggplot(tdat, aes(Sex, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-sex.pdf"), p, height = 2, width = 6)
p

Figure 2: Mortality by sex.

Oxygen

Relationship oxygen to outcome

tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    supp_oxy2 = factor(supp_oxy2, labels = c("Not required", "Required")), 
    D28_death) %>%
  group_by(supp_oxy2) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(supp_oxy2, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", x = "Supplemental oxygen")
p2 <- ggplot(tdat, aes(supp_oxy2, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", x = "Supplemental oxygen")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-oxygen.pdf"), p, height = 2, width = 6)
p

Figure 3: Mortality by oxygen.

Country

Relationship country to outcome

tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(Country = fct_infreq(PT_CountryName), D28_death) %>%
  group_by(Country) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country.pdf"), p, height = 2, width = 6)
p

Figure 4: Mortality by country.

Site

Relationship site to outcome

tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    D28_death) %>%
  group_by(Country, Site) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()
p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Day 28 mortality by site within country.

Calendar Time

Relationship calendar date to outcome

caldat <- all_dat %>% 
  filter_fas_itt() %>%
  dplyr::count(D28_death, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.15)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Relationship between calendar date and the primary outcome.

Analyses

FAS-ITT

Code

res <- fit_primary_model(fas_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.88	(0.10, 8.22)	1.68 (2.79)	0.55
Intermediate-dose	0.79	(0.39, 1.60)	0.85 (0.32)	0.74
Low-dose with aspirin	0.80	(0.34, 1.82)	0.87 (0.39)	0.70
Therapeutic-dose	3.95	(1.14, 13.02)	4.75 (3.26)	0.02
Ineligible aspirin	4.69	(1.15, 16.80)	5.79 (4.24)	0.02
Age ≥ 60	2.00	(1.05, 3.78)	2.10 (0.70)	0.02
Female	0.34	(0.16, 0.69)	0.37 (0.14)	1.00
Oxygen requirement	3.23	(1.67, 6.25)	3.42 (1.19)	0.00
Australia/New Zealand	0.53	(0.11, 2.36)	0.71 (0.62)	0.80
Nepal	2.51	(0.55, 9.80)	3.17 (2.52)	0.11

Posterior summaries for model parameters (fixed-effects), FAS-ITT.

Odds ratio summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Code

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8726 0.8759 0.8284 0.8378 0.9078 0.8129 0.8534 0.8569

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["gamma_epoch"])

Posterior predictive

Code

y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc))
grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(D28_death),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}
ppc_A <- grp_ppc(sym("AAssignment"))
ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("Country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = Country, grp = site) %>%
  summarise(
    mean_y = mean(D28_death),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )
p0 <- plot_grp_ppc(ppc_A, "Antiviral\nintervention") + labs(x = "")
p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p0 / p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-3-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Posterior predictive distribution versus observed proportion (red diamond).

ACS-ITT

Fit primary model

res <- fit_primary_model(acs_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.83	(0.08, 8.54)	1.68 (3.11)	0.56
Intermediate-dose	0.80	(0.39, 1.64)	0.85 (0.33)	0.73
Low-dose with aspirin	0.82	(0.35, 1.85)	0.89 (0.39)	0.68
Therapeutic-dose	4.00	(1.18, 13.20)	4.83 (3.31)	0.01
Ineligible aspirin	4.48	(1.15, 16.34)	5.57 (4.16)	0.02
Age ≥ 60	2.04	(1.08, 3.85)	2.15 (0.72)	0.01
Female	0.34	(0.16, 0.69)	0.36 (0.14)	1.00
Oxygen requirement	3.24	(1.68, 6.36)	3.44 (1.20)	0.00
Australia/New Zealand	0.54	(0.11, 2.44)	0.72 (0.64)	0.79
Nepal	2.55	(0.54, 10.10)	3.24 (2.65)	0.11

Posterior summaries for model parameters (fixed-effects), ACS-ITT.

Odds ratio summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Code

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8030 0.7854 0.8184 0.8520 0.8222 0.8149 0.7771 0.7611

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["gamma_epoch"])

AVS-ITT

Pre-specified

• excludes epoch

Code

res <- fit_primary_model(
  avs_itt_nona_dat, 
  logisitc_site, 
  ctr = contr.equalprior, 
  outcome = "D28_death", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile", "ctry"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-avs-itt-summary-table-pre-spec")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.96	(0.08, 11.25)	2.15 (4.91)	0.51
Intermediate-dose	1.00	(0.08, 12.46)	2.29 (4.91)	0.50
Low-dose with aspirin	1.22	(0.08, 16.99)	2.99 (6.36)	0.44
Therapeutic-dose	1.15	(0.08, 15.97)	2.85 (6.71)	0.46
Age ≥ 60	0.37	(0.00, 14.12)	2.06 (8.65)	0.69
Female	0.29	(0.00, 9.14)	1.37 (6.11)	0.75
Oxygen requirement	0.15	(0.00, 4.32)	0.64 (2.36)	0.86
CRP (2nd tertile)	0.34	(0.00, 12.50)	1.80 (7.03)	0.71
CRP (3rd tertile)	0.34	(0.00, 12.54)	1.92 (8.69)	0.71
CRP (unknown)	0.49	(0.01, 23.02)	3.38 (19.52)	0.63
Nepal	0.96	(0.14, 6.73)	1.56 (1.89)	0.52

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Code

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9217 0.9981 0.9607 1.0189 0.9518 1.0084 1.0286 1.0321

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["tau_site"])

Reduced Model

• excludes epoch, region, and site

Code

res <- fit_primary_model(
  avs_itt_nona_dat, 
  logistic_mod, 
  ctr = contr.equalprior, 
  outcome = "D28_death", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")

Code

XX <- res$dat$X
XX_m <- colMeans(XX)
XX <- cbind(XX[, 1:6], sweep(XX[, 7:12], 2, XX_m[7:12]))
mdat <- res$dat
mdat$X <- XX
snk <- capture.output(
  mfit <- logistic_mod[["sample"]](
    data = mdat,
    seed = 1234,
    adapt_delta = 0.95,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 4,
    parallel_chains = min(4, parallel::detectCores())
  )
)
mpars <- mfit$metadata()$model_params
keep <- mpars[!grepl("(_raw|epsilon_)", mpars)]
mdrws <- as_draws_rvars(mfit$draws(keep))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$OR <- exp(mdrws$beta[!grepl("(Intercept|rand)", names(mdrws$beta))])
median(mdrws$beta)

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.99	(0.09, 11.38)	2.11 (3.82)	0.50
Intermediate-dose	1.00	(0.08, 12.95)	2.33 (4.53)	0.50
Low-dose with aspirin	1.21	(0.09, 16.23)	2.88 (6.31)	0.45
Therapeutic-dose	1.16	(0.09, 15.85)	2.83 (5.98)	0.46
Age ≥ 60	0.38	(0.00, 14.50)	2.13 (8.86)	0.68
Female	0.30	(0.00, 8.97)	1.47 (9.25)	0.74
Oxygen requirement	0.16	(0.00, 4.59)	0.71 (4.17)	0.85
CRP (2nd tertile)	0.35	(0.00, 13.59)	1.95 (7.31)	0.71
CRP (3rd tertile)	0.34	(0.00, 12.07)	1.79 (7.00)	0.71
CRP (unknown)	0.50	(0.01, 22.88)	3.26 (16.63)	0.63

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Code

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0087 0.9952 1.0113 1.0508 0.9884 1.0704 1.0109 1.0048

Trace plots

Code

mcmc_trace(res$drws["beta"])

End of script.